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ABSTRACT 

The discovery of s hort-period Neptu ne-mass objects, now including the remark- 
able system HD69830 lLovis et al.l (|2006h with three Neptune analogues, raises dif- 
ficult questions about current formation models which may require a global treat- 
ment of the protoplanetary disc. Several formation scenarios have been proposed, 
where most com bine the canonical oligarchic picture of core accretion with type 
I mi gration (e.g. iTerquem fc Papaloizoul l2007t) and planetary atmosphere physics 
(e.g. lAlibert et alj [20061 ). To date, due in part to the computational challenges in- 
volved, published studies have considered only a very small number of progenitors at 
late times. This leaves unaddressed important questions about the global viability of 
the models. We seek to determine whether the most natural model - namely, taking 
the canonical oligarchic picture of core accretion and introducing type I migration - 
can succeed in forming objects of 10 Earth masses and more in the innermost parts 
of the disc. 

This problem is investigated using both traditional semianalytic methods for 
modelling oligarchic growth as well as a new parallel multi-zone N-body code de- 
signed specifically for tre ating planetary formation problems with large dynamic range 
(|McNeil fc Nelsonll2009l ). We find that it is extremely difficult for oligarchic tidal mi- 
gration models to reproduce the observed distribution. Even under many variations of 
the typical parameters, including cases in which after the amount of mass in our disc 
is greatly increased above the standard Hayashi minimum-mass model, we form no 
objects of mass greater than 8 Earth masses. By comparison, it is relatively straight- 
forward to form icy super-Earths. 

We conclude that either the initial conditions of the protoplanetary discs in short- 
period Neptune systems were substantially different from the standard disc models we 
used, or there is important physics yet to be understood and included in models of 
the type we have presented here. 
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1 INTRODUCTION 

At present, there are ~19 known extrasolar planet^] with es- 
timated masses between 0.03 and 0.12 Mj up , or between ~ 10 
and ~40 M®. With one exception - OGLE-05-169L b, at 
2.8 AU - all of the planets have semimajor axes smaller than 
1 AU, and so are reasonably called hot Neptunes. In fact, 
with only one more exception, HD69830 d, all have semima- 
jor axes < 0.23 AU. To determine whether the objects are 
genuine Neptunes, i.e. ice giants, and not merely very large 
rocky bodies better thought of as super-Earths, knowledge 
of the mean densities is required. Three of these objects have 
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known radii, and two have mean de nsities compatible with 
being a Neptune-like body. GJ436 b (iButler et alj|2004 ) has 
a densit y of ~ 1.69 g/cm 3 ^Torres et alJl200Sl h and HAT- 
P-ll b (|Bakos et aiTl2009h has a density of ~ 1.33 g/cm 3 ; 
compare Neptune with 1.64 g/cm 3 . The HD69830 system 
jLovis et all l2006h is particularly interesting. It contains 
three Neptune-like objects: HD69830 b at 0.0785 AU, of 10.5 
M® and eccentricity 0.1; HD69380 c at 0.186 AU, of 12.1 
M® and eccentricity 0.13; and HD69380 d at 0.63 AU, of 
18.4 M®, and eccentricity 0.07 (±0.07, so a near-circular 
orbit is possible). The system was also reported to have a 
disc of warm infrared-emitting dus t located between plan- 
ets c and d (|Beichman et alj l2005h . presumably due to a 
collisionally activ e remnant asteroid belt. More recent work 
jLisse et all 12007m suggests instead that the infrared excess 
is due to a debris disc outside the known planets at ~ 1 AU, 
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probably resulting from the breakup of an asteroid. This 
system clearly provides a rigorous test of planet formation 
and migration theories. 

Indeed, short-period Neptune-mass bodies have several 
advantages as probes of planetary origins. They are large 
enough to be observable, but small enough that we need not 
consider gravitational disc instability as a formation pro- 
cess. They are also large enough to undergo significant type 
I migration, but not so large that they can open a gap in 
the d i sc and undergo typ e II migration |Papaloizou fc Linl 
1 19841 : iBrvden et aTTl 19991 : ICrida et all l2006h . Accordingly, 
hot Neptunes can yield less ambiguous tests of type I mi- 
gration than more massive planets which could have signif- 
icantly perturbed the gas disc they were embedded in, and 
their existence provides strong evidence that some kind of 
type I migration is operating in protoplanetary discs. 

Hot Neptunes are therefore useful for exploring the in- 
tersection of the clas sical oligarchic core accretion picture 
dKokubo fc Idalll998l ; e.g. IChambers|[200ll ; iThommes et all 
I2003T ) with type I migration (Ward 1997). Two main classes 
of scenario exist in the literature: one which concentrates on 
the atmospheric gas physics, and one which concentrates on 
the dynamics of the protoplanetary interactions. 

lAlibert et af] (2006) incorporate sophisticated atmo- 
spheric physics and follow the evolution of effectively iso- 
lated cores through the disc as they grow via planetesimal 
accretion, migrate inwards due to type I effects, accrete gas 
after the accretion rate of solids drops, and finally have their 
atmospheric mass reduced after arriving in the short-period 
region by evaporatio n. However, th e histo ry for the HD69830 
system proposed in lAlibert et al.l (I2006T) is difficult to rec - 
oncile with the oligarchic paradigm ( Kokubo fc Ida|[l998h . 
Their best-fitting model involves exactly three seed objects 
of masses M = 0.6 at semimajor axes 3 AU, 6.5 AU, 
and 8 AU, which migrate to small semimajor axis through a 
mostly pristine planetesim al disc and therefore can accrete 
a fair amount of material l|Tanaka fc Idalll999l ). Where did 
these three progenitor objects come from? In an oligarchic 
framework, accretion i n a narrow region - even in the pres- 
ence of type I drag (see | Kominami et al . 2005l: lMcNeir et al.l 
120051 : iDaisaka et al.ll2006l ) - results in many roughly equal- 
mass bodies separated by a distance of 10 ~ 20 Hill radii. 
These bodies then interact and merge in the giant-impact 
phase of planet formation, and consume the accessible rem- 
nant planetesimal disc. That is, if a half-Earth-mass seed 
can successfully form at 3 AU, there should be a large num- 
ber of similar seeds of varying mass inside, each of which 
will accrete and migrate in similar ways. Moreover, the in- 
terior seeds are likely to have formed first. The net result 
is that any such inward-migrating seed should be migrating 
through a region highly depleted by the previous seeds, and 
by the time an 0.6 Mm core is formed at 3 AU much of the 
interior disc should hav e gone to compl etion, and possibly 
formed its own planets (|Chambers!l2008l ). There is no obvi- 
ous mechanism to suppress embryo growth everywhere in the 
disc except at three specific locations. An additional problem 
is that simulations have demonstrated that cores migrating 
through a planetes imal swarm are unab le to grow at the rate 
prescribed by the lAlibert et alj (|2006l ) model. In the pres- 
ence of gas drag, mean motion resonances cause the majority 
of the interior planetesimals to be shepherded rather than 
being accreted, resulting in planets whose masses are too 



low, and in the wrong ratio, compared t o those observed in 
the HD69830 system (|Pavne et alj|2009h . 

Though not aim i ng a t HD69830 in particular, 
iTerquem fc Papaloizoul (120071 ) study the formation of hot 
super-Earths and Neptunes by following the evolution of 
10 — 25 planets of 0.1 or 1 Mj, placed interior to 2 AU under 
type I drag, include tidal interactions with the star, and use 
an inner cavity in the gas disc at ~ 0.05 AU. For various 
disc parameters, they succeed in making several objects of 
mass ^ 8M9, with a maximum of 12 Mfg. They find that 
the typical result has between 2-5 planets, usually on near- 
commensurable orbits (with strict commensurability often 
being broken by tidal circularization) . It is not clear whether 
this will scale up to larger masses, as they performed only 
one run with total mass larger than 12 Mq (namely 25 Mq), 
and so an HD69830-like system would not appear in their 
results even if the model would have succeeded in producing 
them. Very little material was lost from their runs, suggest- 
ing it might be possible. However, their initial configurations 
are difficult to reconcile with an oligarchic migration process. 

Scenarios involving oligarchic formation and type I mi- 
gration do not exhaust the possible formation histories of 
the low-mass hot exoplanet population, alt h ough they are 
arguably the most natural. I Raymond et al.l (|2008l ) surveys 
several other proposed possibilities (for the terrestrial-mass 
regime): in situ formation; shepherding by migrating giant 
planets or secular resonance sweeping; tidal circularization; 
and photoevaporation of giant planets. We will concentrate 
instead on the simplest oligarchic type I migration picture, 
which should have more success self-consistently generating 
an icy planet population, and attempt to determine whether 
the fiducial models can reproduce the observed distribution 
of hot Neptunes. If they can, all the better; if they cannot, 
then the specifics of their failure may point in the direc- 
tion of a solution and help us choose between the various 
possibilities on offer. 

To understand how short-period Neptunes are formed, 
we need to move toward self-consistent global N-body mod- 
els such as those which have proved useful in understanding 
the formation of the terrestrial planets and the outer solar 
system. One problem presents itself immediately: the short- 
period exoplanet problem has a formation time to dynamical 
time ratio ~ 30 times larger than the equivalent terrestrial 
formation problem, and ~ 1000 times that of the equivalent 
Jovian core formation problem. The lifetime of the gas disc 
(with a probable upper limit of ~ 6 Myr) is a common ref- 
erence time-scale for each problem, as gas giants must form 
while there is still gas for them to accrete, and even ice gi- 
ants which attain short-period orbits must have migrated 
while there is gas present. However, the characteristic or- 
bital periods for each formation problem vary from 0.01 yr 
at 0.05 AU to 0.35 yr at 0.5 AU to 11 yr at 5.0 AU. (In- 
deed, 0.01 yr is optimistic, as many exoplanets are closer to 
their parent stars than 0.05 AU.) Since most state-of-the-art 
planetary N-body codes require the integration timestep to 
be some small fraction of the orbital period of the inner- 
most object, preserving the same wall-clock run times that 
researchers have become accustomed to would usually re- 
quire reducing the number of planetesimals in a simulation 
by orders of magnitude, producing an unacceptable decrease 
in resolution. 

It is extremely unlikely that this difficulty can be elim- 
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inated entirely, as it is a consequence of the strong depen- 
dence of orbital period on semimajor axis in the Kepler prob- 
lem. That said, the challenge can be managed to some ex- 
tent by taking advantage of the scale separation caused by 
the troublesome dynamic range. The standard approaches 
use a common drift timestep for all particles (and therefore 
"over-integrate" the more distant objects) and compute the 
(non-encountering) forces between the particles at the same 
frequency, both of which involve far more computation on 
the distant, slow- moving objects than is necessary to pre- 
serve qualitatively accurate dynamics. In I McNeil fc Nelson! 
(2009 ) the authors comb i ne various techniques in the litera- 
ture IjDuncan et al ]| 19981 IChamber s 1999; S aha fc Tremaind 
1 19921 ) to construct a new algorithm which allows for radial 
zones with different timesteps and different inter-zone force 
ev aluation freque n cies, b ut re duces to th e prove n techniques 
of [5 uncan et all (|l998i ") and IChambers! ((l999l ) for objects 
within the same zone. This allows new trade-offs between 

force accuracy and run ti me. 

In our first paper, iMcNeil fc Nelson! (2009), we ad- 
dressed the numerical challenges of studying oligarchic mod- 
els of short-period exoplanet formation. Here we apply a new 
code with a parallel implementation of those methods, to de- 
termine the "reference population" of planets predicted by 
the fiducial models. We integrate global planetesimal discs 
extending from 0.05 AU to 10 AU under various models of 
the protoplanetary disc. Our models, both semianalytic and 
N-body, are described in ^2j and the simulation conditions 
in |3] Simulation results are presented in |4] and discussed 
in SjS] and we conclude in SJS] 



2 DESCRIPTION OF MODELS 

To model the formation of planets in a global disc with a 
large dynamic range, we will divide the system into three 
stages: (1) the first stage, corresponding to the first 0.4 
Myr after our nominal starting time t=0, which we will 
model using a s e mian alytic treatment based upon that of 
iThommes et al.1 (120031 ); (2) the middle stage, from t=0.4 
Myr to 6 Myr, which we will model using the new multiscale 
N-body code; and (3) the late post-gas stage, from 6 Myr 
to 100 Myr, during which we will evolve the inner regions of 
the disc using a traditional SyMBA implementation. 

As usual, the semianalytic model serves two purposes. 
It allows for rough exploration of parameter space and pro- 
vides self-consistent initial conditions for the second stage. 
In practice, we find that for the purposes of serving as initial 
conditions, the model is often more than is necessary: as long 
as the initial mass of the growing protoplanets is much less 
than their final mass, and they have many encounters before 
they undergo significant migration, the systems do not ap- 
pear to show strong dependence on the details of the initial 
conditions. As a quasi-equilibrium, oligarchy is quite robust. 
Nevertheless, although using the model does not eliminate 
the problems caused by initial conditions with inconsistent 
histories, it does help mitigate them. 

To explain our approach, we will first introduce the gen- 
eral features of the gaseous and solid material of the proto- 
planetary disc in sections 12.11 and 12.21 and then describe 
the prescriptions for the aerodynamic and type I drags ap- 
plied to objects by the gas disc in sections l2.3l and [2. 41 These 



models are shared by both our semianalytic approach for the 
early stage, summarized in !2.5l and our N-body approach to 
the later stages, summarized in[ 



2.1 Gas disc 

For simplicity we will consider only models resembli ng those 
of the minimum-mass solar nebula (MMSN) of lHavashil 
l|l98ll ). in which all physical disc parameters may be ex- 
pressed as simple functions of the cylindrical radius from 
the star r and the height above the disc midplane z. 
We take the volume density of the gas to be 

p gas (r, z) = Egas/v^zo exp(- z 2 /2zl ) (1) 

where (r, z) are cylindrical coordinates and zo is the disc 
thickness. We set 



Sf AU (r/AU)" 



(2) 



where £f AU is the gas surface density at 1 AU (1704 g/cm 2 
in the MMSN) and a gives the radial dependence of the 
density (a = 1.5 in the MMSN). In practice, the quantity of 
interest is usually the ratio of mass in the disc to the mass 
in the MMSN. We label this ratio / en ii, and determine the 
appropriate E® AU by normalizing so that the amount of gas 
mass from 0.05 AU to 15 AU in each simulation is / en h times 
the MMSN value. In general terms the disc masses that we 
adopt are consist ent with disc masses inferred fro m sub- 
mm observ ations |Andrews fc Williams] 120051 . 120071 ). More 
specifically, lAndrews fc Williams! (|2005h indicate that more 
than one third of discs in the Taurus- Auriga region have 
masses which exceed that of the MMSN model, with simi- 
lar statis tics applying to discs observ ed in the p Ophiuchus 
complex ([Andrews fc Williams! 120071 ). 
For the disc height we take 



zo 



0.0472 AU (r/AU) 5/4 



(3) 



Note that this power-law in r is traditional in the N-body 
planetary formation community but differs from the canon- 
ical choice in the planetary hydrodynamics community of 
taking a constant zo/r ratio of 0.05 or 0.07. (Although the 
difference may appear minor, it changes the surface density 
power-law at which neighbouring equal-mass objects under- 
going type I migration switch from convergent to divergent 
migration, and is known to be important when reconciling 
results fro m different simulations inv olving capture into res- 
onance: see lCresswell fc N elson 20o3.) We assume the gas is 
in perfect cylindrical rotation. We parametrize the pressure 
support by 



77 = 0.6(z /ffl) 
with typical r\ ~ 0.001, and take 



W 1 ~ 2 V 



(4) 



(•») 



where «k op is the orbital velocity for a circular orbit at r. 

The dissipation of the gas disc with time t is treated 
by introducing a single uniform exponential damping time- 
scale, 

p gas oc exp(-t/r dccay ) (6) 

following iKominami fc Ida! (|2002h . We neglect the gravita- 
tional potential due to the disc. 
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2.2 Solid material 



2.4 Type I migration 



At time t — 0, we set a fixed gas-to-rock ratio of E gas = 
240 E roc k (from the classical values of ~ 7 a/cm 2 for solids 
at 1 AU and ~1700g/cm 2 for the gas. lHavashlll98lf ), and 
therefore take the initial surface density in rocky material 
to be of the form 



,(r/AU)- 



(7) 



where r is the cylindrical radi us and EVatt is the su rface 
density of solids at 1 AU. After iThommes et ail (120031 ). we 
introduce a smoothed snow line beyond which the amount 
of material in solids is enhanced due to the temperature de- 
creasing sufficiently to allow the condensation of ices. Plac- 
ing the snow line at £?i oc and using a smoothing scale for 
the transition of S S m = 0.25 AU, and using an enhancement 
factor Sonh, we write 

S = (0.5 tanh((r - Si oc )/S sm ) + 0.5) (8) 

so that the total initial surface density in solids 

S so lid = E roc k(l — S) -+- 5 C nhS roc k(S) (9) 

The form of S is unimportant, and is chosen simply to soften 
the discontinuity for numerical purposes. 

We assume that all solid material has a density of 2.0 
g/cm 3 . 



2.3 Aerodynamic drag 

We will apply an aerodynamic drag to all bodies (although 
the effects will usually be negligible on objects larger than 
1000 km). We take the drag time-scale as 



8 Pm r m 1 

3 pgas Cd Vrcl 



(10) 



where p m is the mass density of an object, r m its radius, 
w re l its velocity relative to the gas velocity at its position, 
and for simplicity we use a drag efficiency of Cd = 1. In our 
semianalytic model, we will use the approximately equiva- 
lent form 



1 



e m (Cd/2) Pgas a Q 



(11) 



where m is the planetesimal's mass, a its semimajor axis, 
e m its eccentricity, and £1 its orbital frequency. 

Orbit-averaging the above expression to determine the 
planctesimal migr ation rate v m in th e small e, i (inclina- 
tion), and n limi t, lAdachi et alj l|l976t ) (after correction by 
iKarv et al.lfl993l ) 



da I 

dt I aero 



find 
-2 



Tacr- 



v + 



5 2 . 1 -2 , 2 
o Z 

a 



1/2 



_5_ 

4 + 16 



2 _i_ _ -2 



(12) 



For the N-body code, we will use the acceleration 

v - v sas 



a a 



dv 
~dt 



(13) 



For se mianalytic purposes, w e use the type I migration equa- 
tion of ITanaka et alj <|2002l ). which gives an orbit-averaged 
migration rate dr/dt of 



VM 



-Co(2.7+l.la) 



M 



M. 



© 



M, 



© 



rO. (14) 



with Taero from eq. 1101 



for an object of mass M at distance r with orbital frequency 
Q, where we introduce c a as a parameter to incorporate un- 
certainty in the migration efficiency. We note that significant 
contributions to the corotation torque may arise due to non- 
linea r effects, resulting i n a significantly reduced migratio n 
rate (|Masset et alj|2006t IPaardekooper fc Papaloizou]|2009l ) . 
such that c a < 1. The semianalytic model will assume that 
embryos are always on circular orbits and so we do not need 
an expression for the eccentricity damping. 

For N-body calculations, we will use the full i nstan - 
tane ous specific force expre ssions of ITanaka et alj (|2002l ) 
and pTanaka fc Ward! (2004), as given in appendix A of 
iDaisaka et al.1 (|2006l ) . which we do not repeat here. Our only 
modification is to introduce a c a parameter into the relevant 
radial migration term of the form of eq. 1141 

There have been a number of developments over recent 
years which have led to modifications of the basic picture of 
how type I migration may operate under differing assump- 
tions about the underlying protoplanetary disc structure. In 
addition to the corotation torque effects mentioned above, 
it has been noted that regions in the disc where the surface 
density profile has a positive gradient (such that the surface 
density increases outward) may act as planet traps, where 
planet ary migration may not occur at all (jMasset et al.l 
l2006bl ). A planet trap may exist near to the star where the 
stellar magneto sphere clears a low density cavity in the disc 
|Lin et al.l ll996). As discussed in the introduction to this pa- 
per, and in section [5^21 other researchers have considered the 
role of this cavity in the formation of short-period Neptune 
and super-Earth planets, and not surprisingly have found 
that the accumulation of planetary embryos there can lead 
to the formation of planets of the required mass. However, 
the formation of multiple Neptune and/or super- Earth-mass 
planet systems with a broad range of semimajor axes, such 
as HD69830, do not arise naturally from these models. For- 
mation of a multiple planet system inside the cavity is possi- 
ble, but the gravitational potential well of the star likely pre- 
vents scattering out to orbital radii that wou ld be required 
to explain HD69830. iMorbidelli et al.1 (|2008h examined the 
possibility of planetary growth at a planet trap located near 
the snowline. Probably the most plausible explanation for 
a planet trap being located there is the model of inside- 
out disc cl earance by magnetohydrodynam ic turbulence sug- 
ge sted by IChiang fc Murrav-Clavl l|2007h . The simulations 
by IMorbidelli et al.l (|2008l ) indicate that multiple planets at 
a planet trap located at the snowline can indeed undergo 
close encounters and collisions through mutual scattering, 
allowing large bodies to grow there. In addition, scattering 
out of the trap was observed, allowing planets to migrate 
inward after growth through tidal interaction with the in- 
terior disc. One uncertainty which remains in this picture 
is how clear of material the inner disc needs to be f or the 
inside-out clearing model of IChiang fc Murray- Clavl (|2007h 
to operate. X-rays generated in the stellar corona need to be 
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able to penetrate deep into the disc and out to large radii 
for the planet trap to be located at the snowline, leading to 
uncertainty about how much type I migration can actually 
ensue once planets are scattered out of the trap. Given this 
uncertainty in the model, we have chosen to focus on the 
most generic type I migration scenario in the present study, 
which does not include planet traps located at specific radii. 



2.5 Semianalytic approach 

We follow iMcNeil et alJ (|2005h . which is derived from 
ommes et al.l l|2003h . replacing only the formula for type 
I drag used there (due to lPapaloizou fc Larwoodll200C[ ) with 
that of Tanaka and Ward feq. I14p . We will forego repeating 
the derivations and simply describe the resulting equations. 

We simulate the oligarchic migration formation scenario 
using a continuous two-component model consisting of pro- 
toplanetary embryos of mass M(a) which accrete mass from 
a planetesimal field distribution with surface density S m (a) 
orbiting a star of mass Mq, with gravitational constant 
G. The embryos (of density pm) are assumed to be kept 
at a constant fixed separation b in single-planet Hill units 
(rH = (M/3M0) 1//3 a) as a consequence of the usual oli- 
garchic equilibration between the increased separation due 
to scattering and the decreased separati on due to accretion 
i|Kokubo fc Idalll998T l. (As discussed in lMcNeil et~aT1l2005l . 
this approximation is questionable at late times when strong 
migration is operating or the discs are very massive, but is 
reasonable during the early stages.) The embryo eccentrici- 
ties are neglected as they are likely to be much smaller than 
those of the planetesimals due to both dynamical friction 
and type I damping (itself a kind of dynamical friction with 
the gas). The planetesimals (of density p m ) are assumed to 
have a uniform mass m (and radius r m ), and at a given 
semimajor axis all have an equilibrium eccentricity found 
by equating the time-scale for stirr ing by the embry os by 
damping by aerodynamic drag. See IChambersl ([2006) for a 
detailed description of the weaknesses of the equilibrium ec- 
centricity assumption, and also note that by neglecting the 
contribution of embryo-embryo mergers we underestimate 
the accretion rate. 

The physics in the model is simple, although the alge- 
bra is somewhat tedious. Embryos accrete mass from the 
planetesimals according to the expressiorjf] 



dM \ 
dt I 



3.93 M /6 G 1 ' 2 E m M 2 ' 3 6 2 / 5 C 2 J 5 p 2 g L! 



(15) 



p 1 / 3 oVid mVis ptl 15 

which corresponds to a decrease in planetesimal surface den 
sity 



dt 



-M Q /3 



dM 



3 2 / 3 bira 2 M 1 / 3 dt 



(16) 



Planetesimals migrate due to aerodynamic drag with a ra- 
dial rate v m given by eq. 1121 and embryos migrate due to 
type I effects with a radial rate vm given by eq. 1141 These 
four equations - 1151 1161 1121 and [14] - are numerically inte- 
grated. 



2 The factor of 6 2 / 5 is missing from eq. 4 of lMcNeil et al.l l|2005t ): 
the error was typographical. 



We will assume that the separation b = 10 and that 
all planetesimals have r m = 10 km. We will take the initial 
seed mass for the embryos as M(a,t = 0) = 1.5 x 1O~ 3 M0, 
which has the effect of artificially accelerating the growth of 
the more distant embryos during the earliest stages. 

2.6 N-body approach 

The semianalytic Eulerian approach of the previous section 
can be useful as a rough estimate of the behaviour of the 
system, but contains very limited information about the 
dynamics involved. For later stages in which the interac- 
tions between the embryos are significant a particle-based 
Lagrangian approach using an N-body code is necessary. 

Currently, the most robust N-body algorithm for oli- 
garchic simulations of planet formation on long time- 
scales remains SyMBA {5 uncan et al.lll998h . which derives 
from the original mixed-v ariab le symplectic integrato rs of 
IWisdom fc Holmanl (|l99ll ) and iKinoshita et all (l99ll ) but 
has been improved to treat close encounters between par- 
ticles. Unfortunately, it requires a common base timestep 
for all particles, which is set to some small fraction (typi- 
cally ~ 1/15 or less) of the inner most period so that peri- 
centre passage is always resolved. |Levison fc Duncan! l2OO0l 
introduced a variant which could handle occasional objects 
crossing the usual innermost boundary by smoothly switch- 
ing to a Bulirsch-Stoer integrator, following up on the in- 
novations of Chambers 1999, but it becomes impractical 
whe n boundary crossings ar e common.) In a companion pa- 
per ([McNeil fc Nelson! l2009l ) we introduce a new algorithm 
NAOKO which allows for multiple radial zones with distinct 
timesteps and can vary the number of force evaluations be- 
tween different zones, making possible a new trade-off be- 
tween the force accuracy between distant objects and speed. 
We have implemented a parallel version of NAOKO in the 
planet formation code MIRANDA, which is basically a par- 
allel SyMBA implementation. 

As in the semianalytics, we have two classes of objects, 
embryos and planetesimals, where the embryos can merge 
with each other and with the planetesimals, but the plan- 
etesimals do not self-interact either gravitationally or col- 
lisionally. Two objects (assumed spherical) merge if their 
physical separation is less than the sum of their physical 
radii, and form one new object, conserving mass and lin- 
ear momentum. No fra gmentation is consid ered. Following 
standard practice (e.g. lThommes et al.ll2003h we use "super- 
planetesimals" , and replace the planetesimals of mass m 
we would prefer to use with larger objects of mass m sp 
which behave as objects of mass m with respect to non- 
gravitational interactions like gas drag, and serve as rep- 
resentative elements of an underlying planetesimal popula- 
tion. As long as the number of super-planetesimals is large 
and the embryo mass is greater than the super-planetesimal 
mass (M S> m sp ), the accretion rate is insensitive to this 
approximation. (Taking a uniform m for the planetesimals 
and neglecting their self-interactions are considerably more 
damaging simplifications than using super-planetesimals in 
any event.) 

Interactions with the gas disc are limited to aerody- 
namic drag and type I drag, as described in sections 12.31 
and !2.4l The super-planetesimals will experience an aerody- 
namic drag corresponding not to the physical radius of the 
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integrated tracer, but to that of the underlying planetesimal 
(typically ~ 10 km). 

For the computationally challenging second stage dur- 
ing which the gas is present, we will integr ate the system us- 
ing M IRANDA in its new NAOKO mode l|McNeil fc Nelson! 
2009). This allows the disc to be divided into distinct radial 
zones, and objects in each zone are integrated using differ- 
ent time steps. In the simulations presented here we used 
four zones, with timesteps chosen so that all objects had at 
least ~ 15 steps per orbit: namely, At = 0.53 yr for objects 
outside 4 AU; At = 0.13 yr for 1.6-4 AU; At = 0.013 yr for 
0.34-1.6 AU; and At = 0.00083 yr for 0.05-0.34 AU. This cor- 
responded to timestep ratios between outer and inner zones 
of 4, 10, and 16, respectively. These ratios also describe the 
ratios of the frequencies on which the interzone forces were 
evaluated; intrazone forces were evaluated at each (zone) 
step. Boundaries between zones had associated transition 
zones centred there, of widths 0.5 AU, 0.1 AU, and 0.04 AU 
from outermost to innermost, in which the objects smoothly 
experienced both timesteps to avoid artificial kicks in veloc- 
ity. For the late stage after the gas disc has dissipated, the 
number of particles in the inner zone of interest will have de- 
creased enough that we can return to the traditional SyMBA 
method to simulate the final giant-impact stage, with a fixed 
timestep of 0.0007 yr. 



3 SIMULATION CONDITIONS 

Computational power being limited, it was not feasible even 
with the new code to perform multiple runs for each param- 
eter set of interest. Instead, we made compromises between 
coverage of parameter space and reproducibility of each run, 
and between concentrating on physically plausible scenarios 
and less plausible but informative limiting cases. Table [T] 
lists the resulting choices, and where two simulation labels 
are given we ran two instantiations which differed only in 
the random number seeds used to define the initial Keple- 
rian angular variables of the particles. 

We used mass enhancement factors of / cnn = 3, 5, 10, 
to cover the range from the likely enhancements above the 
MMSN needed to form the solar system to something con- 
siderably larger. (Again, £f Au is determined from normaliz- 
ing the amount of mass between 0.05 AU and 15.0 AU to be 
/ cn h times the MMSN value.) The surface density power-law 
was chosen from a = 1.0, 0.5, 0.001, so that the discs are all 
much flatter than the MMSN (with a = 1.5). The MMSN 
value was skipped in the production runs presented here be- 
cause preliminary low resolution simulations indicated that 
disc models with a high degree of central mass concentra- 
tion do not successfully form surviving short-period plan- 
etary systems with significant amounts of mass. This may 
be due in part to the fact that such a disc model induces 
divergent type I migration for neighbouring bodies of equal 
mass, as discussed below, thereby reducing the planetary 
growth rate. In add it ion, v iscous disc models based on the 
IShakura fc Sunvaevl (Il973l ) 'alpha' prescription for viscos- 
ity tend to generate shallower surface density distri butions 
jPapaloizou fc Terquemlll999l ; iFogg fc Nelsonl 120071 ). Note 
that with zo/r oc r 1//4 (eq. |3J), the type I dr/dt rate is inde- 
pendent of a for equal M at a = 1.0, and shows convergent 
migration for a < 1.0. We used disc dissipation time-scales 



Table 1. Simulation parameters. 



Simulation 


/cnh 


Sioc [AU] 


a 


Ca 


Tdccay [Myr] 


S01A, S01B 


3 


2.7 


0.001 


0.3 


1 


S02A 


3 


2.7 


0.001 


0.3 


2 


S03A, S03B 


3 


2.7 


0.001 


1.0 


1 


S04A 


3 


2.7 


0.001 


1.0 


2 


S05A, S05B 


3 


2.7 


0.5 


0.3 


1 


S06A 


3 


2.7 


0.5 


0.3 


2 


S07A, S07B 


3 


2.7 


0.5 


1.0 


1 


S08A 


3 


2.7 


0.5 


1.0 


2 


S09A, S09B 


3 


2.7 


1.0 


0.3 


1 


S10A, S10B 


3 


2.7 


1.0 


1.0 


1 


SUA 


3 


2.7 


1.0 


1.0 


2 


S12A, S12B 


5 


2.7 


0.001 


0.3 


1 


S13A 


5 


2.7 


0.001 


0.3 


2 


SUA, S14B 


5 


2.7 


0.001 


1.0 


1 


S15A 


5 


2.7 


0.001 


1.0 


2 


S16A, S16B 


5 


2.7 


0.5 


0.3 


1 


S17A, S17B 


5 


2.7 


0.5 


1.0 


1 


S18A 


5 


2.7 


0.5 


1.0 


2 


S19A, S19B 


5 


2.7 


1.0 


1.0 


1 


S20A 


10 


2.7 


0.5 


0.3 


1 


S21A 


10 


2.7 


0.5 


1.0 


1 


S22A 


10 


2.7 


1.0 


0.3 


1 


S23A 


3 


4 


0.001 


0.3 


1 


S24A 


3 


4 


0.001 


1.0 


1 


S25A 


3 


4 


0.5 


0.3 


1 


S26A 


3 


4 


0.5 


1.0 


1 


S27A 


3 


4 


1.0 


0.3 


1 


S28A 


3 


4 


1.0 


1.0 


1 


S29A 


5 


4 


0.001 


0.3 


1 


S30A 


5 


4 


0.001 


1.0 


1 


S31A 


5 


4 


0.5 


0.3 


1 


S32A 


5 


4 


0.5 


1.0 


1 


S33A 


5 


4 


1.0 


0.3 


1 


S34A 


5 


4 


1.0 


1.0 


1 


S35A 


10 


4 


0.5 


1.0 


1 


S36A 


10 


4 


1.0 


0.3 


1 


S37A 


10 


4 


1.0 


1.0 


1 



(eq. [6]) of 1 Myr and 2 Myr, which are roughly compatible 
with the o bserved disc decay times inferred from observa- 
tions fe.g. lHaisch et al.|[200lh . The migration efficiency c a 
was set to either 0.3 or 1.0. The snow line was placed at ei- 
ther Sioc — 2.7 AU or 4.0 AU, with an enhancement factor of 
S C nh = 4. Some st udies suggest that a snow line closer to the 
star (1.6-1.8 AU, iLecar et all 120061 ) and a lower enhance- 



ment factor (~ 2.2, Loddersl " 20031) may be more realistic 



but others argue that a m uch higher enhancement can oc- 
cur IjCiesla fc C uzzi 2006). Since many of the Neptune- mass 
objects found to date orbit lower-luminosity stars, a closer- 
in snow line may more accurately model the environments 
of the known short-period Neptune planets. However, to fa- 
cilitate comparison with other work we will keep the mass 
of the central star at Mq and use the historical snow line. 
(The S\oc = 4.0 AU runs were an experiment motivated by 
some preliminary simulations which suggested that delay- 
ing the onset of embryo growth past the snow line could 
hel p prevent promi sing objects from falling into the star; 
c.f. IChambersfl2006l .) 

Each simulation was evolved for the first 0.4 Myr using 
the semianalytic model of H2.5I This yielded a distribution of 
embryo mass M(a) and planetesimal surface density E m (a). 
These were then discretized into an N-body particle distribu- 
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tion extending from 1 AU to ~ 10 AU, with spacing between 
the embryos fixed at b ~ 10 and the super-planetesimal 
mass m sp chosen to ensure that locally M/m ap ^ 5 with 
r m = 10 km . (This is a somewh at more relaxed condition 
than used in iMcNeil et all (|2005l ). which required the rela- 
tionship to hold globally and used one uniform m sp .) The 
resulting particle sets had between 40 and 64 fully inter- 
acting embryos, with > 10000 super-planetesimals for the 
main S\ oc = 2.7 AU runs and ~4000 super-planetesimals for 
the lower-resolution Si oc = 4.0 AU experiments. The initial 
embryo eccentricities and inclinations were arbitrarily set 
to 0.001 and 0.0005, respectively, and the planetesimal ec- 
centricities to their (semianalytically-estimated) equilibrium 
values. The simulations rapidly reach their random-velocity 
equilibrium. The particle distributions were then integrated 
until t=6 Myr using the NAOKO algorithm of US] This 6 
Myr time-scale was chosen so that even in the simulations 
with Tdecay = 2 Myr, at least 95% of the original gas mass 
is gone. Finally, the resulting protoplanets inside of 2 AU 
were then integrated until t=100 Myr using the traditional 
SyMBA algorithm. (That is, the remnant planetesimal disc 
was removed entirely, as were all embryos with semimajor 
axis > 2 AU.) 

Each integration took roughly 3 ~ 4 weeks of runtime 
on an 8-processor node for the 0.4-6 Myr phase, and then 
another week to two weeks in serial mode for the 6-100 Myr 
phase. Without the use of both parallelism and the new 
NAOKO algorithm, the simulations would have taken an 
impractically long time. 



4 RESULTS 

In this section we present some sample runs in H4.ll includ- 
ing both representative cases showing the radial evolution 
of the embryos in H4.1.1I and the evolution of a few of our 
"successful" runs in H4.1.2l and ij4.1.3l We discuss the global 
results in H4.2I giving an overview in H4.2.1I The planetary 
mass distributions are covered in H4.2.2I the distribution of 
ices in H4.2.3I various statistics on the resulting configura- 
tions in £14.31 the final number of planets and the radial mix- 
ing in i]4.3.1l the total surviving mass in the inner region in 
^ 14.3.21 and, finally, the possibility of future mergers and the 
long-term stability of the systems in i]4.3.3l 



4.1 Description of individual runs 

4-1.1 Sample migration histories 

Figurcs[T][2] and|3]offer example overviews of the radial evo- 
lution of the embryos. In S14A, a flat disc (a = 0.001) with 
moderate enhancement / cn h = 5, the outer regions are still 
chaotic at 6 Myr. In S07A, a steeper disc with a — 0.5 and 
low enhancement / cn h = 3, an object from beyond the snow 
line grows massive enough sufficiently early to migrate in- 
wards much faster than its neighbours and compress the ear- 
lier objects into a much smaller radial region, causing both 
mergers and occasional outward ejections (e.g. just after 2 
Myr). There is little communication between the innermost 
and the outermost regions, creating a large empty region in 
(a,t) space devoid of planetary embryos. In S19A, a — 1.0 
with moderate enhancement / cn h = 5, we see that almost all 



Figure 1. Semimajor axis, perihelion, and aphelion for each em- 
bryo over time in simulation S14A. 
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of the inner material is rapidly lost to the inner edge of the 
simulation, and much of the material which eventually forms 
the resulting inner system comes from far beyond the snow 
line. However, some of the resulting planets incorporate ma- 
terial which started inside of the snow line and was scattered 
out beyond it before migrating back in. Each simulation had 
c a = 1.0 and Tdecay = 1 Myr, and therefore the majority 
of their tidal migration is complete by 3 Myr, although as 
S14A shows, a considerable amount of radial migration can 
continue due to embryo-embryo interactions long after the 
gas di sc is gone. This can be compared with IMcNeil et al.l 
(2005h 's study of terrestrial planet oligarchic migration sce- 
narios, which predicted a tripartite division into an interior 
region of strong convoying behavior (where planets rapidly 
lock into mean-motion resonance and migrate in tandem), 
a transition region where objects 'slide' towards their final 
destinations, and an outer region which remains chaotic. In 
their study, they looked only at the ct = 1.5 and a = 1.0 
cases, and did not consider any a generating convergent mi- 
gration. 



4.1.2 Run S09A 

The run S09A produced the single most massive object of 
any simulation, a planet of 7.2 Mm. This particular run 
had the parameters / cn h = 3, a = 1.0, c a = 0.3, Td ccay = 
1 Myr, and Si oc = 2.7 AU. At t=0.4 Myr, the simulation 
start, there were roughly comparable amounts of mass in 
embryos and planetesimals inside of 2 AU, M to t = l-l M0and 
m to t = 1.4 Mm, respectively. By 1 Myr, both have increased, 
to M to t=2.3 M0and m tc .t = 1.8 VU. The amount of embryo 
material in the region gradually increases to 11.1 M0at 6 
Myr. This increase is non-monotonic, as embryos are occa- 
sionally lost in two ways. First, some are scattered beyond 
the inner simulation edge at 0.05 AU where they are re- 
moved. Second, some are scattered beyond the 0.05-2 AU 
window, in which case they will remain in the system but 
no longer contribute to the mass (at least unless type I 
migration or embryo-embryo scattering brings them back 
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Figure 2. Scmimajor axis, perihelion, and aphelion for each em- 
bryo over time in simulation S07A. 



< 




Figure 4. Run S09A planets interior to 2 AU at t=6 Myr, planet 
mass versus semimajor axis, with bars indicating eccentric excur- 
sion. 



Figure 3. Scmimajor axis, perihelion, and aphelion for each em- 
bryo over time in simulation S19A. 




3 4 
time [Myr] 

in again). By contrast, the planetesimal material reaches a 
maximum of ~ 2.3M0 at 2.2 Myr, as it migrates into the 
region via both Jacobi shepherding by embryos and aerody- 
namic drag, and then falls off to ~ 1.5M0 for the remainder 
of the simulation, as material is consumed but new material 
is brought in to replace it. 

As shown in figure |4j at 6 Myr there are nine em- 
bryos with a < 2 AU, and they fall into two distinct cat- 
egories: four objects inside of 0.8 AU, with masses less 
than 1.2 which vary by factors of three, and five well- 
spaced objects exterior to 0.8 AU which all have masses 
~ 1.6 — 1.7M0, where the inner three have startlingly sim- 
ilar masses, namely 1.608 Mg, 1.609 Mq, and 1.607 M^. 
Figure [5] shows that this is simply a coincidental result of 
equilibration processes, as the objects follow roughly similar 
growth curves - an early period of fast accretion transition- 
ing to a much slower phase after 3 Myr when the gas is 
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Figure 5. Run S09A planet mass versus time for all objects. 
From t=6 Myr, only the embryos interior to 2 AU are followed. 




time [Myr] 



mostly gone - but at times the objects differ in mass by 
a factor of ~ 2. The three ~ 1.61M0 objects also come 
from very different locations in the disc: the outermost em- 
bryo started its life at 2.4 AU, the middle at 3.3 AU, and 
the innermost at 1.6 AU. During the giant impact phase, the 
objects at 0.87 and 1.06 AU merge quite quickly, at ~7 Myr, 
as do the two innermost bodies. By ~ 40 Myr, the original 
(t=6 Myr) 1.7 M0body has consumed three more objects 
to reach its terminal mass, while the remaining objects were 
scattered outside 2 AU. 
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Figure 6. Run S16B planets at t=6 Myr. 



Figure 9. Set of final planets at t = 100 Myr, planet mass versus 
semimajor axis. Pies as in fig. [7] 
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4.1.3 Run S16B 

The run S16B produced the largest amount of surviving 
embryo mass interior to 2 AU at 100 Myr of all the sim- 
ulations, 16.75 Mm. This particular run had the parameters 
/onh = 5, a = 0.5, c a = 0.3, T dccay = 1 Myr, and Si oc = 2.7 
AU. This region started with M to t = l-00 M0and m to t = l-3 
Mq, reaching final values of M to t = 16.8 M0and m to t = l-9 
M e at 6 Myr. 

As shown in figure [6] there are 8 embryos with at least 
an Earth mass of material (including a hot Earth of 1.01 
M0at 0.37 AU). The planet masses vary more than in S09A, 
although it is still true that the planets with the lowest 
masses are on the interior and the highest masses have semi- 
major axes between 1.2 and 1.6 AU. The system is basically 
quiescent until a burst of activity between 14 and 15 Myr, 
during which the planet at 1 AU merges with the planets at 
0.7 AU and 1.9 AU, and the planet at 0.6 AU merges with 
the planets at 0.8 AU and 1.6 AU. After this, the system 
does nothing but exchange angular momentum with little 
variation in semimajor axis. 

Comparing S09A with S16B suggests that maximizing 
the mass of the final planet does not require maximizing the 
mass of the planets when the giant impact phase begins; 
both simulations produced maximum objects of mass 6 — 
7M0, despite S16B having 1.5 times the mass of S09A at 6 
Myr. 



4.2 Synthesis of simulation outcomes 

4-2.1 Overview 

Snapshots of the planetary systems which resulted at T=100 
Myr are shown in figures [7] and [8] with each planet being 
represented by a pie showing its fraction of rock and ice, 
and with horizontal bars indicating the eccentric excursion 
(and not, as is often done, ±5 Hill radii.) It is worth not- 
ing that the maximum ice fraction possible is 0.75. Bodies 
composed entirely of material whose provenance lay interior 
to the snow line are 100% rock. Several general features of 



0.5 1.0 1.5 

semimajor axis [AU] 



the simulation outcomes are immediately apparent in these 
figures. 

Most of the objects, especially the smallest ones, have 
modest eccentricity. S07B and S17B are exceptions, in which 
the outermost body has a nonneglible e (0.16 and 0.29, re- 
spectively.) There is also the case of S24A, which is mani- 
festly unstable. The two outer planets had been exchanging 
a and e since an order-swapping encounter at 70 Myr, and 
had suffered a series of particularly close encounters ~ 5 
Myr before the end of the simulation. It might have been 
expected that all of the larger objects would have the small- 
est eccentricities due to dynamical friction and type I drag, 
but once the gas disc has dissipated and the local planetes- 
imals have been consumed, the remaining large planets can 
generate significant eccentricities. 

It is also clear that stochasticity is playing a very im- 
portant role. While in some cases the same initial conditions 
generate similar final configurations even after the variation 
of the angles (e.g. S01A/B, S05A/B), in other cases the re- 
sulting systems have almost no resemblance to each other 
(e.g. S09A/B, S14A/B, S17A/B). This confirms the regret- 
table fact that the simulation outcomes are sufficiently sen- 
sitive to the details of the encounter and merger history that 
isolated runs convey limited information, although they can 
certainly demonstrate that a given scenario is possible. 

Figure[9]shows the complete collection of resulting plan- 
ets (combined from all simulations). The majority of objects 
we form are small: the median planet mass is 1.07 Mm, and 
90% of the objects have masses below 3.55 Mm. Defining 
success as the production of an object with greater than 4 
M0interior to 2 AU, there were only eight successful runs: 
S04A, S06A, S09A, S13A (2 objects), S16A, S16B (2), S17B 
(2), and S33A. Only one object (in S33A) came from a run 
with Sioc = 4 AU, and so our early, mildly encouraging ex- 
periments with a more distant snow line proved unfruitful 
in the production runs. 5 of the 7 successes (not double- 
counting S16) had /enh = 5, 2 had / cnh = 3, and - sig- 
nificantly - none had / cn h = 10. 5 of the 7 successes had 
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Figure 7. Final planetary systems at t = 100 Myr by simulation, planet mass in Earth masses versus semimajor axis in AU. Each planet 
is plotted as a pie with radius proportional to its physical radius, with the blue (brown) slice corresponding to the fraction of its mass in 
ice (rock). Horizontal bars describe the eccentric excursion, i.e. they extend from a(l — e) to a(l + e). The labels specify the simulation 
reference number as listed in table [l] 



c a — 0.3. Every a value generated at least one success (2, 
3, 2 for a = 0.001,0.5,1.0, respectively), and likewise for 
Tdccay (4 for 1 Myr, 3 for 2 Myr). The single most success- 
ful parameter set was that associated with S16A and S16B, 
namely / cnh = 5, a = 0.5, c a = 0.3, r dccay = 1 Myr, 5 loc = 
2.7 AU. 

It bears noting that the combined plots and histograms 
that we present in this section should be interpreted as 'pop- 
ulation synthesis' only in a loose sense. There have been no 
corrections made for the shape of the parameter space, which 
includes duplications of some parameter runs, and which 
was not chosen to match any expected distribution of disc 
parameters in actual protoplanetary systems. Accordingly, 
they should be taken merely as summaries of our particular 
outcomes, and not as predictions for what any real-world 
observer would see. 



4-2.2 Planetary mass distribution 

Keeping that warning in mind, figure [10] shows the number of 
bodies of a given mass produced across all the simulations. 



Figure 10. Summary of planet mass outcomes. 
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Figure 8. Final planets at t = 100 Myr; as in fig, 171 



As mentioned previously, half of the planets formed have 
mass ^ . Beyond 3 , the numbers fall off dramat- 

ically. This is in agreement with some (unpublished) early 
exploratory low-resolution runs which had great difficulty 
forming planets greater than 4M0 regardless of the disc en- 
hancement used. The median mass reached in the simula- 
tions is a full order of magnitude - 17 times - smaller than 
Neptune. 

We can also consider not only the successful runs (in 
which we have small-number statistics problems) but look 
at the complete dataset: see table [2] Perhaps unsurpris- 
ingly, the greater the enhancement factor, the greater the 
median mass. However, this increase in the median value 
does not correspond to an increase in the maximum; indeed, 
the greater the enhancement / cn h, the smaller the resulting 
maximum mass of a planet. A more robust measure is the 
change in the mean mass of the top quintile, which shows 
that / cn h = 5 produced more large bodies than / onn = 10. 
A similar wrap-around occurs with a, where a = 0.5 has a 
higher median and top-quintile-mean than either a — 0.001 
or a — 1.0. The existence of such maxima in parameter 
space suggests that some natural methods for increasing the 
planet growth rate - such as simply increasing the initial 
surface density of the disc - may not help increase the final 



mass, as there are feedback mechanisms operating which re- 
sist the formation of larger bodies (namely the rapid inward 
migration of massive bodies which form early in a gas-rich 
environment). In part, this result could merely be statisti- 
cal noise, as of the 7 simulations with / cnn = 10, only 3 
produced planets (and all of those had weak migration with 
c a — 0.3), but this low number is itself an e xample of the 
proble m. This agrees with the predictions of iDaisaka et al.l 
and (except at v ery high disc masses) with the re- 
sults of lChamberj <|2008h . 



Table [2] reveals several other (apparent) correlations be- 
tween the simulation parameters and the resulting planet 
masses beyond those involving / cnn and a. Decreasing the 
migration efficiency (c a ) improved the final mass, changing 
the median mass from I.OMq) to 1.4M0 and the top quin- 
tile mean from 3.OM0 to 4.5M0. The use of a more distant 
Sloe = 4, contrary to some of our early experiments, tended 
to reduce the planet mass. Finally, the runs with Tdccay = 
2 Myr had smaller masses than those with Td eca y = 1 Myr, 
suggesting that the ability to bring in mass from farther out 
in the disc can overcome the loss of material to the star (or, 
at least, the simulation edge.) However, given the wide vari- 
ation in outcomes for the same parameter set - and our use 
of no more than two instantiations of each - these results 
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Table 2. Summary of resulting planets by parameter; values are 
taken over all planets produced in the subset of runs with param- 



eter = value. 


Parameter 


Value 




planetary mass 


statistics [Mm] 






Median Maximum 


Mean of top quintile 


/cull 


3 


0.826 


7.194 


2.667 


/enh 


5 


1.455 


6.359 


4.428 


/enh 


10 


2.058 


3.917 


3.917 


a 


0.001 


1.025 


5.000 


3.353 


a 


0.5 


1.412 


6.359 


4.478 


a 


1.0 


0.987 


7.194 


2.878 


"Sloe 


2.7 AU 


1.124 


7.194 


3.980 


"Sloe 


4.0 AU 


0.982 


5.193 


2.901 


Ca 


0.3 


1.430 


7.194 


4.482 


Ca 


1.0 


1.014 


6.359 


3.030 


^dccay 


1 Myr 


1.044 


7.194 


3.447 


^decay 


2 Myr 


1.250 


5.000 


4.402 



should be treated cautiously, despite our averaging over the 
suite. 

In some runs, a considerable amount of mass was re- 
moved from the simulation by falling off our inner boundary 
of 0.05 AU. In the / cnh =10, Si oc =2.7 AU runs S20A, S21A, 
and S22A, the lost planets totalled 34 M0(with the maxi- 
mum mass of a lost planet being 5.2 Mfg), 38 Mjj (maximum 
4.0 M(g), and 67 Mjj (maximum 7.1 Mq). (The continuous, 
fixed-spacing semianalytics generating the initial conditions 
are probably unreliable at such a large enhancement.) In the 
/enh = 3 and / cnn = 5 cases, the amount of mass lost var- 
ied from none (in all the discs with a = 0.001 except S15A, 
which had c a = 1 and T dccay = 2 Myr) to 21-24 M^in S18A 
and S19A/B. Excluding the / cn h=10 cases, none of the lost 
bodies were larger than 3.3 M0(in S16A), and the median 
mass of a body which was lost was 1.9 Mq. 

4-2.3 Ice fractions 

Figure ^demonstrates that larger-mass objects tend to have 
higher ice fractions, and given our initial assumption that 
75% of the material past the snow line was in ices, the ice 
fraction is a proxy for the radial transport. The median ice 
fraction for objects greater than 1 is 0.60, and for ob- 
jects smaller, 0.36. There were only 7 objects which con- 
tained no ices at all (1 from S01A; 3 from S23A; 2 from 
S24A; 2 from S28A), and all but one of those had masses 
smaller than 0.5 Mm, the exception having 0.75 Mm. Ev- 
ery ice-free body was formed from a run with / enn = 3 and 
Tdccay = 1 Myr, suggesting that we should only expect to 
find completely rocky bodies in low-enhancement discs with 
short disc lifetimes. Moreover, all but one ice-free planet was 
from a simulation with S\ oc = 4 AU, which suggests that in 
the more physically realistic scenarios we should expect al- 
most no bodies which fail to contain material from beyond 
the snow line. (Of course, 'ice-free' and 'completely rocky' 
here are relative to the assignment of ice fractions at t=0.4 
Myr, when our N-body integrations began, although they 
should be assigned in a way consistent with the snow en- 
hancement.) It is notable, as shown in figure 1111 that the 
entire range of possible ice fractions is covered, from bod- 
ies which are completely ice-free to bodies which reach the 
maximal 75% value. 

Figure[l2]shows that the lower the ice fraction, the more 



Figure 11. Summary of ice fraction outcomes. 




fraction of material in ice 



of an object's mass was consumed (or started) as an embryo, 
and not a planetesimal. At the start of the simulation, all 
embryos are either inside the snow line, in which case we as- 
sume they are composed entirely of rock, or outside the snow 
line, in which case we assume they are composed of 75% ice 
and 25% rock. (We neglect the early semianalytic evolution 
of the models before t = 0.4 Myr; as a result, we slightly 
underestimate both the initial ice fractions of embryos in 
the inner regions and their consumption of planetesimals.) 
Accordingly, at t=0.4 Myr, all embryos are either located at 
(1, 0) or (1, 0.75) on this plot: they are composed entirely of 
embryo material, and have ice fractions of either or 0.75. 

At first, objects located in both regions accrete local 
material from both embryos and planetesimals, which de- 
creases their embryo fractions, but leaves their ice fractions 
constant as they accrete from nearby material which has 
the same ice fractions as they do (whether or 0.75). Both 
interior and exterior objects therefore move left on the dia- 
gram. Eventually migration becomes important, which tends 
to move ice-rich material into the inner regions. Accretion 
then raises the ice fractions of the rocky bodies (which now 
have icy material to consume) and lowers the ice fractions 
of the icy bodies (which previously had mostly icy material 
available but now find more rocky material in their feeding 
zone). Since so much material is brought into the 0.05-2 AU 
zone from outside, this tends to concentrate the resulting 
planets at high ice fraction with roughly 20-40% of their 
mass coming from embryos (the median embryo mass frac- 
tion is 0.28). Almost all objects with very low embryo mass 
fraction (< 0.2), which therefore consumed almost all of 
their material in planetesimals, have very high ice fraction. 
These values may be biased somewhat by our choice of in- 
ner edge, which artificially depletes the innermost regions of 
rocky planetesimal material. 

In the innermost regions, a < 0.25 AU, somewhat 
counter-intuitively, we find smaller bodies with higher ice 
content. All objects save one (an 0.43 MQice-free object 
from S28A) had masses 1 Mjj < M < 4 Mq , and were com- 
posed of at least 50% ice. Excluding the rocky outlier, the 
median ice fraction was 0.69. Since the maximum ice frac- 
tion possible in our simulations is 0.75, these objects are 
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Figure 12. Ice fraction versus amount of material consumed by 
an object in the form of an embryo (as opposed to a planetesimal). 
By construction all objects have ice fractions SC 0.75. 



5 0.4 




0.2 0.4 0.6 0.8 

fraction of material sourced as embryo 

composed almost entirely of material from beyond the snow 
line. 



4.3 Simulation statistics 

IChambersl (|200ll ) (§4, which we follow closely here) intro- 
duced a useful set of dimensionless statistics for comparing 
the results of terrestrial planet simulations, but as he notes 
they are of general applicability. (Given our inability to build 
planets larger than 8 Mm, the statistics are rather more 
closely applicable than we had hoped.) His set includes: 

(i) N, the number of bodies (here planets). 

(ii) Sm, the fraction of the total mass in the largest object. 

(iii) S s , a spacing statist ic, loosely related to t he empir- 
ical instability time-scal e l|Chambers et al.1 1 19961 ; see also 
llwasaki fc Ohtsukil [200rj for a more recent investigation of 
the problem, which suggests that the dependence on the 
mass should be closer to 0.29 ~ 2/7), defined by 

1/4 



N - 



+ a, 



where a m i n and a max are the limiting semimajor axes (here 
set to 0.05 AU and 2.0 AU), and m is the mean planet mass. 

(iv) Sd, the normalized angular momentum deficit, de- 
fined by 



Sd = 



ft 1 " ^/l-e|cosij) 



(18) 



with summa tion over the planets, indexed by j. (See, e.g.. 
lLaskarlll997h 

(v) S c , a concentration statistic, given by 



S c = max 



Em,- 



Em, [\og 10 (a/aj)Y 



(19) 



where we take the maximum value the argument reaches 
over all values of the variable a in the interval [0.05 AU, 
2 AU]. The higher the value, the more 'concentrated' the 



Table 3. Summary statistics for each run; definitions as in i|4,3l 
For comparison, the label "TERR" describes the four terrestrial 
planets in our own system, and 3NEP the three Neptunes in the 
HD69830 system. Undefined values are labelled N/A, and un- 
known values left blank. Simulations S21A, S35A, and S37A re- 
sulted in no planets interior to 2 AU, and so are not listed. 



Sim. 


N 


Mtot 


TTl 


Sm 


S s 


s d 


S c 


S r 


S01A 


3 


1.50 


0.50 


0.55 


90.23 


0.008 


73.55 


0.44 


S01B 


4 


1.51 


0.38 


0.31 


64.54 


0.003 


79. 20 


0.46 


q° 2 "a 




3.12 


1.56 


0.59 


135.79 


0.014 


51.42 


1 . 3 5 


b03A 


: 

5 


2. Of) 


0.41 


0.33 


47.35 


0.OO2 


8 4.02 


(-'.'> 9 


S 1. 1 - 1 13 


- 


1.90 


0.95 


0.80 


153.74 


0.018 


5 3 . i 8 


. 9 • 


S04A 


1 


4.06 


4.06 


1 .00 


N/A 


0.016 


N/A 


2.54 


S 5 A 


2 


3. 78 


1.89 


0.54 


129.43 


°'°"f 


35.04 


1.49 




3 


4.20 


1.40 


0.39 


'^'l' 'J' 


O.OO^i 


43.67 


1 . ; j 1 




3 


9.61 


' ! ' J " 


0.51 


,)(). / 1 


0.007 


48.89 


3.0.) 


a "A 
?° - , 


6 


5.17 


0.80 


0.10 


31.49 


2.1c-04 


"!'!'-- 


2.99 


bO t B 


- 


5.02 


2.81 


0.03 


117.20 


0.014 


33. ,-r> 


2.31 


SOS A 


8 


5. :j6 


0.0 9 


0.22 


23. i 4 


4.5c-05 


7.67 


10.02 


S 9 A 


1 


7.19 


7.19 


1 .00 


N/A 


0.02 3 


N/A 


3.49 


S 9 B 


8 


9.78 


1.22 


0.23 


20.01 


1.4c-04 


20.13 


2.10 


SI OA 


7 


7.84 


1.12 


0. 30 


2 1.58 


l.lc-04 


4. 10 


lO.Oo 


S10B 


6 








31.57 


2.6c-05 




3.07 


SUA 




97 


o 'i8 


I! 






544 68 




S12A 


5 


4.98 


1.00 


0.40 


37.97 


3.8e-04 


84.87 


1.94 


S12B 


4 


4.42 


1. 10 


0.36 


49.3 3 


0.001 


75.57 


1 .83 


S13A 


5 


14.92 


2.98 


0.34 


28.86 


9.4e-05 


23.58 


5.62 


S14A 


2 


5.92 


2.96 


0.60 


115.67 


0.007 


88.91 


2.82 


S14B 


6 


9.67 


1.61 


0.20 


26.93 


6.2e-04 


36.31 


2.73 


S15A 


6 


15.16 


2.53 


0.25 


24.07 


0.007 


6.86 


16.94 


S16A 


4 


10.29 


2.57 


0.42 


39.93 


2.7e-04 


54.90 


2.53 


S16B 


4 


16.75 


4.19 


0.37 


35.35 


0.012 


21.20 


3.39 


S17A 


5 


8.47 


1.09 


0.44 


33.25 


2.1e-04 


2.21 


30.17 


S17B 


2 


11.36 


5.68 


0.56 


98.28 


0.041 


13.34 


6.12 


S18A 


2 


1.01 


0.50 


0.67 


180.19 


5.0e-07 


32.16 


3.74 


S19A 


6 


6.02 


1.00 


0.21 


30.32 


1.0e-06 


6.10 


14.58 


S19B 


2 


2.05 


1.03 


0.52 


150.72 


9.3e-06 


1426.51 


11.62 


S20A 


3 


5.42 


1.81 


0.38 


65.44 


2.7e-06 


25.43 


33.28 


S22A 


1 


3.33 


3.33 


1.00 


N/A 


2.4e-06 


N/A 


2.05 


S23A 


3 


1.12 


0.37 


0.67 


96.94 


0.002 


122.19 


0.23 


S24A 


3 


1.27 


0.42 


0.41 


93.96 


0.034 


42.33 


0.45 


S25A 


3 


2.30 


0.77 


0.42 


81.07 


0.002 


94.80 


0.67 


S26A 


2 


2.80 


1.40 


0.64 


139.49 


0.011 


70.28 


1.32 


S27A 


5 


4.89 


0.98 


0.37 


38.14 


0.001 


50.19 


1.28 


S28A 


9 


7.42 


0.82 


0.15 


19.90 


2.1e-04 


11.39 


2.70 


S29A 


3 


3.42 


1.14 


0.43 


73.39 


6.3e-04 


119.09 


1.36 


S30A 


3 


5.42 


1.81 


0.65 


65.42 


9.7e-04 


79.77 


1.89 


S31A 


5 


10.65 


2.13 


0.35 


31.40 


1.4e-04 


47.78 


2.35 


S32A 


4 


7.00 


1.75 


0.34 


43.97 


0.002 


9.02 


7.42 


S33A 


3 


7.85 


2.62 


0.66 


59.65 


0.007 


45.37 


3.39 


S34A 


5 


4.16 


0.83 


0.38 


39.71 


2.2e-04 


19.50 


6.39 


S36A 


1 


3.92 


3.92 


1.00 


N/A 


6.6c-05 


N/A 


2.29 


TERR 4 
3NBP 3 


1.98 
41.00 


0.49 
13.67 


0.51 
0.45 


6.29 
4.12 


0.002 
0.004 


89.49 
7.05 





radial distribution is, although this statistic is not a proxy 
for the mass surface density. 

(vi) S r , a radial mixing statistic, given by 



Sr — 



Em. 



mj\a,j, 



a i.f 



(20) 



where and / refer to the original (at t=0.4 Myr) and final 
(at t=100 Myr) semimajor axes of each object which ulti- 
mately becomes part of o ne of the resultin g planets. This 
statistic was introduced in lChambers! (|200ll ) to quantify the 
degree of mixing induced by gravitational scattering. Being 
linear in the change in semimajor axis, this statistic may 
not prove as useful measuring this effect in simulations with 
high migration rates, but will quantify the degree of type I 
migration instead. 

These statistics, along with the total mass M to t and 
mean mass m, are listed in table [3] as well as the corre- 
sponding values for our terrestrial system (Mercury, Venus, 
Earth, and Mars) and the system of three Neptune analogues 
orbiting HD69830. 



4-3.1 Final number of planets and radial mixing 

A somewhat surprising outcome of the simulations is that 
there were so many systems with high N which survived: 9 
of the 48 simulations resulted in systems with N > 5 even 
after 100 Myr. Some of this surprise is explained by the fact 
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that intuitions trained on the terrestrial regime can fail for 
a < 0.5 AU, where the same radial separation corresponds 
to a much greater Hill separation, and many of the high- 
N systems have multiple objects in the interior. Removing 
these objects from consideration reduces the number of high- 
N systems to 5. That said, these systems (both the original 
9 and the reduced 5) cover every a, include / on h = 3 and 
5, both Tdecay values, both She values, and all but one have 
c a — 1.0. This ubiquity suggests that the unexpectedly large 
number of planets is likely to be the result of strong migra- 
tion selecting stable configurations. 

When applied to our runs, the radial mixing parameter 
S r defined by equation [20] provides a measure of the degree 
of migration experienced by all components of the final plan- 
etary system. For example, a (somewhat unrealistic) system 
of two equal mass planets which begin life with semimajor 
axes a ~ 2 AU, and subsequently migrate inward without 
further accretion to the radial distance ~ 0.2 AU will have 
a value S r — 10 (this value becomes S r — 20 if the final 
stopping location decreases to 0.1 AU; it becomes S r — 1 if 
the final stopping distance increases to 1 AU). As figure [TBI 
shows, while there are many low-N runs with S r comparable 
to the high-N values, there are no high-N runs with low S r . 
Of the 9 N > 5 runs, all had radial mixing S r > 2, with 
4 having S r > 10 (44%); of the 40 N < 5 runs, 40% had 
S r < 2, and only 3 had S r > 10 (7%). This suggests - al- 
beit weakly, given the small numbers involved - that high-N 
systems consist of bodies that formed over a wide range of 
radii and migrated into the interior region. There are also 
no high-N runs with Mtot < 5 Mq (numerous planets add 
up to a significant amount of total mass). The absence of 
any high-N systems with / en h = 10 is likely the result of the 
would-be planets migrating out of our integration region, 
and another example of how in this regime higher surface 
densities can hurt more than they help. Indeed, all three 
runs which resulted in no planets interior to 2 AU (S21A, 
S35A, and S37A) all had / on h = 10 and c a — 1, although 
here our choice of outer boundary may be playing a role; we 
return to this issue in H5.ll 

It is probable that the perturbations from giant planets 
would significantly lower these numbers. Their formation is 
not modelled here: accretion of gas is not treated, n or are 
any objects of a > 2 AU after 6 Mvr. lChambersI l|200ll ) notes 
that one explanation consistent with the observed decrease 
in terrestrial planet number between his simulations and 
those of some previous authors was his incorp oration of the 
Jovian s, and a similar effect should occur here. iMcNeil et al.l 
<|2005h also had to invoke an external random velocity source 
to reduce the number of planets in late-stage prototerrestrial 
systems produced via oligarchic migration. Of course, even 
had we incorporated a gas accretion model, we do not form 
any objects large e nough to serve as useful seeds for the 
iPollack et al.l (1996) picture while the gas is present. 

4-3.2 Total surviving mass 

Of primary importance for our purposes is the total amount 
of mass which survives interior to 2 AU at t=100 Myr. Only 
6 runs succeeded in producing more than 10 M0of material 
inside 2 AU. Despite the small number, it is striking that 
every one had / cn h = 5 and S\ oc — 2.7 AU, with a < 1.0. S16 



Figure 13. Number of planets versus radial mixing parameter. 
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14. Number of planets versus total and mean planet 
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successful, with both instantiations producing over 10 Mg, 
and S16B producing two of the largest bodies in the suite 
and the single largest total mass, 16.75 Mq. Recall that for 
our vertical profile (eq.[3]), a — 1.0 is the transition between 
convergent and divergent migration for equal-mass bodies, 
and that for a < 1.0, neighbouring equal-mass migrating 
objects converge. At the outset of this project we had hoped 
that the compression might increase the resulting masses, 
and both the behaviour at the high Mtot limit and summing 
over the other parameters as in table [2] support this hope. 
When the Mtot limit is lowered to 5 Mgj, 7 of the 23 runs 
(30%) had a — 1.0, and so any advantage that convergent 
values of a possess seems limited to the extremes (which 
is where one might expect to first observe increases in the 
variance of the mass with a lowered a.) 

Figure [14] shows the variation of Mtot and ra with 
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the number of planets. The cluster of high-mass runs with 



M, 



> 



15 Mmhad nearly the median number of planets 



(4 ^ N ^ 6), while both lower and higher values of N had 
lower maximum Mtot (11.4 M0and 9.8 M^g, respectively). 
This hints that median N values produce the largest masses, 
but there are only three high-mass runs involved. Moreover, 
the four runs with N > 6 also have mean and median M to t 
greater than or comparable to the median runs, even if the 
maximum is greatly reduced (unsurprisingly, having more 
planets tends to translate into a greater total mass). It does 
appear that both the maximum value and the variance de- 
crease with small and large N; there are no runs with N ^ 6 
and Mtot < 5 . There is also a strong correlation between 
the number of resulting planets and their mean mass, such 
that the larger N is, the lower m is. The curves which bound 
m as JV varies are quite regular, which given the volume of 
parameter space that our search covered is strong evidence 
that we must make large modifications to our initial condi- 
tions or our physics to increase m to Neptune-like masses at 
N > 30 

One way to see the difficulty is to compare the total 
mass in a run with the mean mass (fig. 1151 zero-mass out- 
comes are suppressed.). There is a cluster of very low-mass 
runs with m ~O.5M0, but at larger m and M to t the region 
of achieved values opens up. Clearly m ^ M to t, which sets 
the lower boundary curve. Removing the high-m outliers, 
the / C nh = 3 and / cn h = 5 cases show similar behaviour, 
except that the / cn h = 5 set reaches higher total masses 
and has a larger variance. Not only do the / en h = 10 cases 
fail to improve on the / cn h = 5 case, each is inferior to 
many / cn h = 3 runs. Worse yet, three zero-mass runs were 
not plotted, although some of those had full migration and 
lost implausible amounts of mass to the star. This leads us 
to tentatively conclude that any mechanism that increases 

3 It is darkly amusing to note that naively extrapolating the best 
fit line for the maxima predicts that the desired m ~ ISMgjoccurs 
at TV ~ —11, a situation which should defeat even the admirable 
ingenuity of today's observational exoplanetary community. 



Figure 17. Mass and separation of neighbouring pairs at t=100 
Myr. 




20 40 60 80 100 

semimajor axis separation in mass-weighted Hill units 

the rate of planetary growth during early times is likely to 
reduce, and not increase, the mass of the largest surviving 
bodies in the presence of significant inward type I migration. 

Figure [T6l shows the fraction of the total mass contained 
in the largest body (by definition, this cannot be smaller 
than 1/N). Slightly under half of the planet-producing sim- 
ulations have S m > 0.5, and the trend toward lower S m with 
larger N is clear. For comparison, the current HD69830 value 
has N = 3, S m = 0.45, well within the achieved range. 

4-3.3 Future mergers and long-term stability 

It is unlikely that most of the resulting systems will un- 
dergo many more mergers (save S24A). Figure [T7] shows 
the mass in each pair of neighbouring planets in all the re- 
sulting systems against the mass-weighted Hill separation 
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Figure 18. Evolution of median separation of planets interior to 
2 AU. 




(6 = ((Mi + M 2 )/3M ) 1/3 (A/ 1 a 1 + M 2 a 2 )/{M 1 + M 2 )) of 
the pair. There are only 5 pairs of planets involving total 
masses greater than 8M^, and 4 of those have b > 20. The 
median b is ~ 22, and the median mass of objects closer 
than the median separation is only 2.3M0 . If we follow the 
mass-separation relationship for each pair through the sim- 
ulations, we find that the vast majority of pairs which have 
small separations have combined masses < 4 (which 
helps explain why we do not form very many large objects) . 

Figure [TS] shows the evolution of the median mass- 
weighted Hill separation of all objects inside 2 AU for all 
simulations until 6 Myr. While the data are quite noisy, it 
is clear that almost all simulations show a definite evolu- 
tion from their initial value of ~ 7 (~ 9 in single-planet 
units) to ~ 15 (~ 19); the median outcome is 16.5. 33 of 
the 48 simulations (69%) have median b > 15, and 14 have 
median b > 20. Only 6 of the runs (12.5%) had median 
b < 13. This difference can be significant as the interaction 
time-scale is a strong fun ction of the separation. Work by 
llwasaki fc Ohtsukil (|2006l ) shows that the time-scale for a 
collection of cold equal-mass objects to suffer an instability 
in the absence of nebular gas when separated by 15-20 mu- 
tual Hill radii can be > 10 10 years. These systems have much 
shorter interaction times in practice due to t he variance 
of th e spacing and the nonzero eccentricities l|Zhou et alj 
120071 ). but the likelihood of the set of planets merging to 
form Neptune-sized objects even when the total mass makes 
it possible is very low. In any case, these total masses are far 
beneath the ~ 41M0 present in the HD69830 system, and 
the maximum mass of 7.19 Mgjis far below the estimated 
18.4 M ffi of HD69830 d. 



5 DISCUSSION 

5.1 Weaknesses of the model 

As usual, many of the approximations invoked to make the 
simulations tractable may have affected the results. 



The semianalytic model used to advance the early 
stages of the model is quite primitive. It uses a fixed Hill 
spa cing of 10 between th e embryos, using neither the origi- 
nal |Kokubc^£jda| (|l998l ) dependence of b on M nor a more 
realistic slower growth. This should be a tolerable error dur- 
ing the first 0.5 Myr for all but the highest-enhancement 
cases. The approximation of the embryo distribution by a 
smooth function means that most of the interembryo dy- 
namics is lost, which can be important when tidal convoys of 
migrating objects locked into mean-motion resonance form. 
In a noisy N-body migration problem, objects which become 
larger than their neighbours due to a merger can push inte- 
rior embryos with it at migration rates larger tha n the naive 
model would predict (as in S07A - see figure [2)l. lChambersl 
(2008) discusses this issue in the context of building semi- 
analytic models of oligarchy. The use of a large seed mass 
Mo makes the discs easier to instantiate in the N-body code 
by lowering the gradient in embryo mass, but at the cost 
of artificially accelerating the growth of the most distant 
objects, which should lead to higher migration. Quite sub- 
stantial growth above this initial mass is required, however, 
to enable trans-snowline embryos to migrate interior to 2 
AU, so this is unlikely to be a significant problem in the 
simulations. 

Enforcing a sharp embryo/planetesimal distinction such 
that planetesimals neither self-gravitate nor accrete, the 
usual procedure in the field among those using Kepler- 
problem symplectic integrators, results in a poor treatment 
of the mass spectrum. No planetesimals can promote them- 
selves to embryos, even if a planetesimal ring forms in which 
the largest body should experience runaway growth and be- 
come a new oligarch. This will suppress embryo growth in 
some cases (when an oligarch should have formed), and in- 
crease it in others (by providing a fresh supply of mate- 
rial f or protoplanets migr ating into the regio n to consume, 
as in iTanaka fc Idal Il999l and the models of lAlibert et al.l 
2006). The use of a uniform characteristic planetesimal ra- 
dius means that we are insensitive to the differing effects 
of aerodynamic drag on different-sized objects. While the 
planetesimal non-accretion problem is difficult to overcome, 
the uniform-mass problem could be handled by varying the 
radius assigned to super-planetesimals, so that some objects 
would trace the behaviour of 1 km bodies, some 10 km, and 
so on. However, in the absence of any collisional physics - 
here, we use a simple hit-and-stick model with no fragmen- 
tation - the small-radius planetesimals will not be replen- 
ished as they would be in a real system. The absence of 
a collisional cloud may increase the eccentricity of the re- 
sulting planets as an object cannot be damped by its own 
debris, but this is probably negligible during the gas-rich 
phase when type I drag is strong and the embryos are on 
effectively circular orbits the majority of the time. 

Ultimately, it will be necessary to use a hy- 
brid N-body/stati st ical scheme along the lines of 
iBromlev fc Kenvonl (|2006l ). and we are currently ex- 
ploring t he incorpo r ation of a variant of the semianalytic 
model of IChambersI (120081 ) into the code. 

The choice of inner and outer edge for the initial solid 
material distribution in the N-body runs - 1.0 AU and ~ 10 
AU - may have changed the outcomes in two ways. In the 
runs with low migration, some objects which should have 
been present in the interior will be missing. However, this 
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is seldom a significant amount of mass. The only disc in 
which there was more than 4M0of material interior to 1 
AU which was excluded was the / en h = 10, a = 1.0 case 
(~8Mgj), which has very high migration; S22A, S36A, and 
S37A all fail to contain any objects inside of 1.5 AU be- 
cause they have all migrated away, and so very little ma- 
terial (if any) would have survived. In the discs with the 
lowest migration rates, where / cn h = 3, only 2.3M0 of mass 
was missed when a — 1.0. (The flatter profiles have more 
of their mass at larger r and so have even less initial mass 
inside of 1 AU.) Therefore the effects on the final results are 
limited: low-enhancement discs which have migration rates 
small enough that we should have included the innermost 
embryos also have too little mass in the region to be in- 
teresting, and high-enhancement discs which have enough 
mass to significantly affect the mass of a final planet have 
migration rates high enough that the missing material is 
long gone by the time the gas disc has dissipated. The inner 
cut-off edge of 0.05 AU will also be responsible for removing 
some particles with small semimajor axis and high eccen- 
tricity, but there are very few of them. On the outer edge, it 
is quite possible that a wider particle disc than 10 AU would 
have resulted in more planets in the high-enhancement runs 
and prevented the zero-planet outcomes (all of which were 
very high-migration configurations with / cn h = 10 and c a = 
1). However, such high enhancements and migration rates 
are at the limits of plausibility to start with, and the initial 
conditions from the semianalytic model are questionable at 
high /cnh- Most importantly, there is no reason to expect 
that the planets which would have been formed would be 
any larger than the ones already formed in the / cn h = 10 
runs, as objects which migrate into the inner regions in such 
runs are often either well-separated or resonantly locked and 
therefore protected against most encounters. 

In our N-body model, the only effect that the gas disc 
has upon the particles is via prescription as a source of drag, 
whether aerodynamic or type I, and we did not include any 
disc physics involving gap opening. This probably had little 
effect, as if we assume a (hydrodynamic-standard) constant 
zo/a ratio of ~ . 05, an d use the gravitothermal condition 
from ICrida et al.1 £2006), gap opening even in an inviscid 
disc requires (3/4)(zo /th) < 1, or th — 3/4H, which gives 
a mass of ~ 53M0 (with "gap opening" defined to be a 10% 
perturbation in surface density). Even assuming a ratio of 
1/2 sufficed to open a gap, we would need ~ 17M0. Things 
become somewhat better in a flared disc (such as the one 
we actually used) with zq/cl oc a 1 ^ 4 , where if we assume 
no kinematic viscosity, an 8 M0 object could possibly start 
opening a gap at ~0.08 AU, near the edge of our simulation 
region. In practice, gap opening is unlikely to be a good 
explanation for the presence of the 'lukewarm' end of the hot 
Neptune population unless it is substantially more efficient 
than currently thought. 

We also neglected the accretion of gas onto the em- 
bryo cores. Given the low planet masses we obtained - 
with median 1 Mm and maximum ~ 8 M(j- the standard 
iPollack et all (|l996h model suggests that even our most mas- 
sive planets would have accreted atmospheres of at most 
1-2 M0. This missing mass should not have affected the 
outcomes significantly, at least directly, as it would have 
changed the Hill radii by less than 10%. However, it is known 
jlnaba fc Ikomall2003l l that growing embyros can have a sig- 



nificantly increased effective capture radius due to the en- 
ergy loss suffered by planetesimals moving through the em- 
bryo atmosphere, which can help overcome the challenge of 
forming giant plan et cores before the gas disc evaporates 
|lnaba et al.l 12003 ). Would including this effect have made 
it possible to form substantially larger cores? 

This question should be considered in a larger context. 
One feature common to many of the new mechanisms being 
discussed, atmosphere-enhanced capture radii included, is 
that they tend to increase the accretion rate at early times. 
Partly this is due to a target bias on the part of scientists 
studying the middle and late stages of planet formation, as 
forming the cores of giant planets on appropriate time-scales 
is an important open problem even in the migration-free 
case, and so efforts have been concentrated in the direc- 
tion of making large objects easier to form. In the absence 
of migration, this is a net positive. The situation is more 
complicated when migration is involved. Increasing the ac- 
cretion rate sounds appealing, but since the type I migration 
rate scales linearly with the mass, it will simultaneously de- 
crease the migration time-scales, and do so at early times 
when the gas density is at its highest. If anything, this will 
make the formation-to-migration time-scale problem worse, 
at leas t in cases where migration is significant. Daisak a et all 
(2006) noted that for a fixed gas-to-dust ratio and gas dis- 
sipation time-scale, increases in initial s urface density can 
result in decreases in final surface density. IChambersl (I2008T ) 
finds that for a disc with alpha viscosity of 0.001, the max- 
imum embryo mass increases as the disc mass increases to 
a high of ~ 3M0 at O.O5M0, and then falls off, so that 
disc masses of O.O3M0 and O.O7M0 both have maximum 
planet mass 0.5 — O.6M0. Including fragmentation does raise 
the maximum mass reached over a wide range of disc mass 
(0.03 - 0.10M Q ), but only to ~ 4M . The case with alpha 
viscosity of 0.01 resembles the case without migration be- 
low disc masses of O.O5M0. Only at very large disc masses, 
greater than 0.10 Mq, are Neptune- like masses recovered. 

Modifying the accretion rate seems unlikely to encour- 
age neighbouring oligarchs to merge at late times, given their 
large separations. These spacings are consequences of the 
high surface density of the discs as well as type I migration. 
Enhanced capture radius mechanisms should have little di- 
rect effect on the giant impact phase, as a body needs to 
encounter roughly its own mass of gaseous material in the 
extended atmosphere to significantly affect its orbit. For ex- 
ample, Earth-mass objects experiencing a near-collision will 
barely notice any atmosphere, however extended. In general, 
changes to the accretion rate in these oligarchic migration 
scenarios should move the window of source material which 
grows and migrates to the innermost regions, but will not 
increase the total mass by the factors of several needed to 
recover a planetary system resembling HD69830 (except pos- 
sibly for massive discs). That is, it should change where the 
planets come from, not how large they are. 

5.2 Comparisons with previous work 

Our results differ significantly from previous results in the 
liter ature, chiefly i n that we fail where they succeed. 

lAlibert et all i|2006l ) succeeded in forming a system 
which resembles the HD69830 system very closely, by re- 
ducing the strength of the migration by a factor of 10 
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and considering only thr ee c ores which migrat e and accrete 
fc.f. lTanaka fcldall 19991 ). As lChambersl (|2008l ) explains, the 
presence of interior embryos with a much shorter dynamical 
time means that inward-migrating cores do not encounter 
a pristine feeding zone, but one which is already substan- 
tially processed. None of our migrating cor e s exp erienced 
growth resembling that in the lAlibert et ail (2006) model. 
IPavne et al.l (|2009h further find that aerodynamic drag can 
significantly decrease the planetesimal accretion rate in such 
scenarios, keeping the migrating core acting as a shepherd 
and not a predator. 

An alternative way of ensuring that a population of hot 
super-Earths or hot Neptunes can form and survive is to in- 
voke an inner disc cavity generated by interaction between 
the disc and the stellar magnetosphere. Such a model has 
been used to explain the existence of hot Jupiters (|Lin et al.l 
19961). and was incorpor ated into the formation models o f 



Brunini fc Cioncol (|2005l ) and lTerquem fc Papaloizoul (|2007l ) 



(who also included tidal inte raction with the centra l star) . 
iBrunini fc Cioncol (|2005l ) and lTerquem fc Papaloizoul (|2007l ) 
both succeeded in forming large short-period Nep tune-mass 
(10-30 M^) cores. Flerauem fc Papaloizoul (|2007l ) form four 
objects of mass ^ 8M 0, the largest being 12M0, and 
IBrunini fc Cion co (2005) form five objects of mass ^ 1OM0, 
with the largest being ~24M0. There are three main dif- 
ferences between their models and ours which likely account 
for the discrepancy. First, bot h are relatively low-resolution: 
iTerquem fc Papaloizoul l|2007l ) use d only 10 Earth-ma s s bod - 
ies in all but one of their runs, and IBrunini fc Cioncol (|2005l ) 
used 100 embryos of 0.5 M0and 200 planetesimals of 0.1 
M0 . Second, and more importantly, it is unlikely that their 
initial masses and spacings can be recovered from a self- 
consistent oligarchic migration model: there is no plausible 
prior configuration which would evolve into (for example) 
12 Earth-mass bodies spaced between 0.1 and 1 AU. Fi- 
nally, and most importantly, in accordance with the cavity 
hypothesis both groups used a disc with an inner edge at 
either 0.05 AU or 0.10 AU, and we did not. Some small 
toy simulations we performed (unpresented here) confirm 
that we can also build bodies of ~ 15M0 by doing so, but 
only relatively close to the boundary. However, from general 
considerations the magnetospheric cavity in this model is 
expected to extend out to a distance of ~ 0.08 AU for an 
assum ed T Tauri star rotation period of 8 days (|Lin et al.l 
1996), causing migration to halt at a distance of ~ 0.05 AU 
from the star. We simply note that there are a number of 
super-Earths and Neptune-mass extrasolar planets orbiting 
with significantly larger semimajor axes (e.g. HD69830 d, or 
OGLE-05-169L b, at 2.8 AU) whose inward migration was 
probably not halted by the presence of an inner disc cavity. 
It is more likely that they, or their progenitor objects, were 
stranded at or near their current location when the gas disc 
dispersed. 

iThommes et alJ <|2007h and I Chambers! l|2008h use semi- 
analytic models, both incorporating atmosphere-enhanced 
capture radii and the latter fragmentati on, and succeed 
in for ming planets of ^ 



lO Mm, although IThommes et al.l 
<|2007h produce far more than IChambersi (|2008l ). They both 
use an approach in which embryos are treated as dis- 
crete but non-interacting objects instead of as a continu- 
ous distribution (rep laci ng the more Eul erian treatment of 
IThommes et al . 2003 and lChambersll2006l with a Lagrangian 



one). In order to maintain the standard oligarchic sepa- 
ration of ~ 10 Hill radii, they merge two embryos when- 
ev er their semimajo r axes differ by less than 7 Hill radii. 
As IChambersi l|2008l) notes, this will result in missing some 
interesting inter-embryo dynamics such as migration con- 
voys (groups of embryos migrat ing in m ean-motion reso- 
nance, as in McNeil et al. 20051 : see also IThommes! liooBI : 
iPapaloizou fc Szuszkiewiczll2005l ). but it will also miss some 
larger and more fundamental dynamics such as deviations 
of the embryo behaviour from the semianalytic predictions, 
and introduce an unphysical orderly merger wave. 

IChambersi l|2008l ) performed a comparison between the 
results of his semianalytic mode l and one of the terre strial- 
planet N-body integrations of iMcNeil et al.l ([2005) , and 
found that the N-body results had both larger median 
masses (O.25M0 vs. O.I6M0) and larger median spacings 
(~ 20 H ill radii vs. 11). Eve n in the original oligarchic 
model of lKokubo fc Idal i|l998l ). there is a weak dependence 
of b on M and E. If in these scenarios the interembryo 
spacing grows to values much larger than 10 while the ob- 
jects stay on n ear-circular orbits, oligar c hic growth can b e 
qu enched 1, 
As lChambersl ( 



Iwas aki fc Ohtsukil 120061 : 
2006) (i 



IZhou et al.1 120071 ). 



notes, in order to maintain a con- 
stant spacing in b within a fixed width, embryos must ac- 
crete 1/3 of their mass via embryo-embryo mergers. This 
accretion mode should vanish if the embryos are too well- 
separated (and too cold, but the models assume the embryos 
are always on circular orbits.) Too large a spacing towards 
the end of the gas phase can also suppress the giant-impact 
phase almost entirely, requiring external sources of stirring 
to form large cores on reasonable time-scales. Accordingly, 
if the interembryo spacing b for large-mass bodies in these 
simulations is closer to 15 — 20 than 10, then calibrating the 
(admittedly merely statistical) prescription for the embryo 
mergers to the no-migration, low-embryo-mass separation 
value of 10 may not be appropriate. 

Despit e the differences, we find some qualitative agree- 
ments with IChambersi ([2008). In his models with migration 
but without fragmentation (the closest match to our runs), 
he finds that the maximum mass of an object is maximized 
at ~ 3M0 for a disc mass of O.O5M0, and that both higher 
and lower disc mass es decrease this number (as we find, and 
iDaisaka et al.1 120061 predicted) . For his simulations which 
include both fragmentation and migration, as mentioned 
above, there is a wide range of disc masses (0.03 — O.IOM0) 
which produce roughly comparable maxima (~ 4M0). The 
only simulations which succeed in getting cores larger than 
1OM0 either fail to include migration or include both frag- 
mentation and migration with a disc mass of > O.IMq, and 
it appears that the only Neptune-like object had a semi- 
major axis outside 2 A U. This is in mild disagreement with 
rTh ommcs et al l ll2007l) who fo und 10 Mmobjects were com- 
mon, which Chambers! (|2008T > argues is the result of their 
choices of a large alpha viscosity, 0.01, and a small planetes- 
imal radius, 1 km, while neglecting fragmentation. 

We suspect that the solution to the short-period Nep- 
tune formation problem lies not in tweaking the accretion 
rate but in modifying the prescription for the interaction 
with the gas disc. Obviously we need some kind of migra- 
tion to build a short-period population, but what properties 
should it have? Ideally, we would prefer a migration which: 
(1) has little to no effect at early times when the gas density 
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is high, so that enhancing the surface density actually re- 
sults in larger embryos and not more embryos lost, (2) when 
active , has a reduced rate from the nominal iTanaka et al.l 
(2002) behaviour, as our c a = 0.3 runs performed notice- 
ably better than our full migration runs, and (3) provides 
some mechanism to encourage large migrating embryos to 
merge instead of locking in convoys or stranding themselves 
at large separations from their neighbours, such as a variable 
migration direction. 

The most promising model of disc-induced migration 
which displays these cha r acter i stics is that presented by 
IPaardekooper fc Mellemal ()2006l . l200sh . in which the migra- 
tion of planets in radiatively inefficient discs was considered. 
This model has the highly desirable property that migration 
is stopped, or even reversed, during early times when the 
disc is optically thick, but inward type I migration is recov- 
ered when the disc density decreases and the gas becomes 
optically thin. An alternative model for the modification of 
type I migration is stochastic migrati on induced by turbu- 
lent density fluctu ations in the disc (jNelson fc Papaloizoul 
12004 : lNelsonll2005l ). Although this does not appear to have 
the same nicely-tailored characteristics for solving the prob- 
lem of short-period Neptune formation as the radiatively 
inefficient migration model, its random walk nature has the 
potential to deliver significant planetary material into the 
interior regions at late times. 



6 CONCLUSIONS 

We performed 48 simulations of various oligarchic migration 
scenarios to determine whether the simplest standard ap- 
proach can succeed in forming a population of short-period 
Neptune systems under common assumptions for the pro- 
toplanetary disc parameters. Multiple numerical techniques 
were applied: semianalytic techniques for the first 0.4 Myr, 
our new parallel multizone N-body code for the accretion 
phase while gas is present up to 6 Myr, and the more tradi- 
tional SyMBA approach for the late stage to f 00 Myr. 

We find that over a wide range of disc conditions, it 
is difficult to form planets of mass greater than 3 — 4Mq. 
Our most successful runs involved ~ 5 times the mass of 



the MMSN, surface density varying as 



-1/2 



a disc decay 



time-scale of 1 Myr, and a migration efficiency of 0.3. Our 
most common planet outcomes are of Earth-mass objects, 
with the terrestrial planets having ice fractions from 0.0 to 
0.75 (the maximum possible in our simulations). The larger 
objects have higher ice fractions, with the median being 
0.60 for objects above 1 Mg, and 0.36 below. In none of 
the cases did we succeed in forming an object of greater 
than 7.5 M0inside 2 AU, much less inside 0.5 AU, and the 
total embryo mass remaining inside 2 AU was always less 
than 17 Mg. The existence of an upper limit and the weak 
dependence o n most parameters i s in accordance with the 
predictions of iDaisaka et al. d2006Tl, and in rough agreement 
with the predictions of Chambers! ! 2008) except at large disc 
masses. Nevertheless, we should be wary of making predic- 
tions based on these results regarding extrasolar planetary 
systems, as they entirely fail to reproduce the short-period 
Neptune planetary population that we know exists. 

Our failure can be compared to several previous suc- 
cesses in the literature, which either (1) adopt initial con- 



ditions which are not easily reconciled with an oligarchic 
growth picture, (2) use an inner edge to the migration (which 
is defensible but will have difficulty explaining more distant 
Neptunes), or (3) neglect inter-embryo dynamics and use an 
embryo merger condition which is calibrated to an effective 
inter-embryo separation (10 Hill radii) which is considerably 
smaller than we observe. 

Varying parameters which we kept constant such as 
the gas-to-dust ratio, incorporating additional accretion 
physics such as fragmentation, moving to extremely large 
disc masses or extremely weak migration, and simply per- 
forming more runs (and hoping for fortuitous late merg- 
ers), could possibly succeed in improving the maximum mass 
reached by a factor of two and therefore into the Neptune- 
like region. It seems quite unlikely that they will increase 
the median mass enough to comfortably produce a popu- 
lation of multiple-planet short-period Neptune systems. We 
conclude that forming a system like HD69830 will probably 
require a significant revision to the simple models explored 
here. 

If the standard oligarchy-plus-type-I-migration picture 
fails to reproduce the observed distribution of short-period 
exoplanets even at more extreme parameter values, then 
we must consider non-standard models. Oligarchy is rela- 
tively well understood both analytically and numerically; 
by comparison type I migration is sensitive to poorly un- 
derstood properties of the gas disc such as disc turbulence 
and local thermodynamic time-scales. In a follow-up pa- 
per we consider the implications for hot exoplanet forma- 
tion via oligarchy of alternate migration models (such as 
IPaardekooper fc Me llcma 2006) which show some promise. 
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